rm(list=ls())
library(MASS)
library(stargazer)
library(knitr)
library(MASS)
library(foreign)
library(lme4)
library(car)
setwd("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\")
fp <- read.csv("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\PIWP_DataSet_2015.csv")
fp <- read.csv("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\PIWP_DataSet_2015copy.csv")
fp <- read.csv("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\PIWP_DataSet_2015 copy.csv")
########Conact with Female Peacekeepers Table 1 Model 1
########  Do you think that female peacekeepers are better than male peacekeepers?
#NOTE: morefemale.unmil.s is the variable without 98 or I dont know
g<-fp$morefemale.unmil.s ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit1 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit1)
coeff.logit <- coefficients(mod.logit1)
provlnp233<-c(1,0,0,1,32.76328,0,0,1,0,0,0,0)
xb.provlnp233<-provlnp233%*%coeff.logit
predict.prob.win <- exp(xb.provlnp233)/(1+exp(xb.provlnp233))
predict.prob.win
# 0.1441589 (0.08402515-0.2381465)
provlnp233<-c(1,1,0,1,32.76328,0,0,1,0,0,0,0)
xb.provlnp233<-provlnp233%*%coeff.logit
predict.prob.win <- exp(xb.provlnp233)/(1+exp(xb.provlnp233))
predict.prob.win
# 0.2680937 (0.15364350-0.4305131)
mod.logit <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit)
mod.logit <- glm(formula=g,data=fp,family=binomial(link="logit"))
Z<-cbind(1,seq(0,1,0.01),0,1,32.76328,0,0,1,0,0,0,0)
Zb <- Z%*%coefficients(mod.logit)
P <- exp(Zb)/(1+exp(Zb))
b.sim <- mvrnorm(1000,coefficients(mod.logit),vcov(mod.logit))
Zb.sim <- Z %*% t(b.sim)
P.sim <- exp(Zb.sim)/(1+exp(Zb.sim))
P.sim.graph <- t(apply(P.sim,1,sort))
P.sim.graph[,975]
P.sim.graph[,25]
########  Do you think that rape is a problem in Liberia?
# rapeproblem.2 is the variable without I dont know or 98
g<-fp$rapeproblem.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado +fp$aflcontact
mod.logit2 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit2)
######## Do you think that women should join the AFL?
######## Do you think that women should join the LNP?
#NOTE:Womenjoinafl.2 and womenjoin.lnp.2 is the variable without  I dont know or 98
g<-fp$Womenjoinafl.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence    + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit3 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit3)
g<-fp$womenjoin.lnp.2 ~ fp$femaleunmilonly + fp$nounmilcontact + fp$Sex +fp$Age + fp$nowvictim + fp$warviolence   + fp$better.2 +fp$armedgroup +fp$femalelnpcontact +fp$Montserrado  +fp$aflcontact
mod.logit4 <- glm(formula=g,data=fp,family=binomial(link="logit"))
summary(mod.logit4)
######## Create Coefficient Table (Table 1 Model 1)
stargazer(mod.logit1,mod.logit2,mod.logit3,mod.logit4, type="html",
dep.var.labels=c("Fem. PK Better","Rape Problem","Women Join Mil.","Women Join Police"),
covariate.labels=c("Contact with Female Peacekeeper","No contact with peacekeeper ","Female","Age","Victim of Crime","Experienced Violence in War","Better off with UNMIL","Armed Group Member", "Contact with Female Police","Born in Monrovia","Contact with AFL") ,
out="models.html", single.row = TRUE,digits = 2, title = "Table 1")
######## Subsets Sample Only: Table 2 and 3
fem <- subset(fp, fp$Sex == "1" , select=c(hcode:whichdialect))
men <- subset(fp, fp$Sex == "0" , select=c(hcode:whichdialect))
with(fem, table(morefemale.unmil, femaleunmilonly))
with(fem, table(morefemale.unmil, maleunmilonly))
with(men, table(morefemale.unmil, femaleunmilonly))
with(men, table(morefemale.unmil, maleunmilonly))
#with(fem, table(morefemale.unmil.s, femaleunmilonly))
#with(fem, table(morefemale.unmil.s, maleunmilonly))
#with(fp, table(morefemale.unmil.s, femaleunmilonly))
#with(fp, table(morefemale.unmil.s, maleunmilonly))
#with(fp, table(Sex, femaleunmilonly))
#with(fp, table(morefemale.unmil, femaleunmilonly))
#with(fp, table(morefemale.unmil, maleunmilonly))
#######  Figure 1: Perceptions of Protection from Rape Provided by Different Actors
#Who do you think can best protect you from being raped?
#rape2fp = prefer female peacekeeper to respond to rape
#rape1mp = prefer male peacekeeper to respond to rape
#rape2femalepolice = prefer female officer to respond to rape
#rape2malepolice = prefer male officer to respond to rape
#rapet2femalecomm = prefer female community leader to respond to rape
#rape2malecomleader = prefer male community leader to respond to rape
summary(as.factor(fp$rape2fp))  #50/1280 = 0.04
summary(as.factor(fp$rape1mp))  #72/1280 = 0.06
summary(as.factor(fp$rape2femalepolice))  #505/1280= 0.39
summary(as.factor(fp$rape2malepolice))  #852/1280 = 0.67
summary(as.factor(fp$rapet2femalecomm)) #259/1280 = 0.20
summary(as.factor(fp$rape2malecomleader)) #534/1280 = 0.42
barplot(matrix(c(0.04,0.06,0.39,0.67,0.20,0.42),nr=2), ylab="Who do you think can best protect you from being raped?â  ", col=c("black","darkgray"),names.arg =c("Peacekeepers","Police","Community Leader"), ylim=c(0,0.8),beside=TRUE)
box()
legend("topleft", c("Female","Male"), pch=15,  col=c("black","darkgray"),  bty="n")
rm(list=ls())
library(MASS)
library(stargazer)
library(knitr)
library(MASS)
library(foreign)
library(lme4)
library(car)
install.packages("knitr")
rm(list=ls())
library(MASS)
library(stargazer)
library(knitr)
library(MASS)
library(foreign)
library(lme4)
library(car)
mon <- read.csv("~/Dropbox/II Solo/Final/Replication/Monrovia2012 (Beber et al).csv")
mon <- read.csv("C:\\Users\\Daniel\\Dropbox\\Journal (II)\\Replication Files\\Karim\\Monrovia2012 (Beber et al).csv")
#wartrauma (q148): Were any years since start of war really hard times for you?
#psu: community fixed effects
summary(as.factor(mon$psu))
summary(as.factor(mon$q19))  #21 per cent had contact in the last year
summary(as.factor(mon$q21)) #  0.848 with a male and  0.152 with female
summary(as.factor(mon$q25))
#Three new columns created using q21 of contact with female peacekeepers (q21_women), contact with male peacekeepers (mon$q21_men), and no contact with peacekeepers (mon$q21_none).
summary(as.factor(mon$q21_women))
summary(as.factor(mon$q21_men))
summary(as.factor(mon$q21_none))
g<-mon$q25 ~ mon$q21_women + mon$q21_none + mon$q1_age +  + mon$sex + mon$savings  + mon$muslim + mon$cognitive + mon$wartrauma + as.factor(mon$psu)
mod.logit1 <- glm(formula=g,data=mon,family=binomial(link="logit"))
summary(mod.logit1)
g<-mon$q25 ~ mon$q21_men + mon$q21_none + mon$q1_age +   mon$sex + mon$savings + mon$muslim + mon$cognitive + mon$wartrauma + as.factor(mon$psu)
mod.logit2 <- glm(formula=g,data=mon,family=binomial(link="logit"))
summary(mod.logit2)
g<-mon$q25 ~  mon$q21_women + mon$q21_men + mon$q1_age  + mon$sex + mon$savings +  mon$muslim + mon$cognitive + mon$wartrauma + as.factor(mon$psu)
mod.logit3 <- glm(formula=g,data=mon,family=binomial(link="logit"))
summary(mod.logit3)
stargazer(mod.logit1, mod.logit2, mod.logit3, type="html",
dep.var.labels=c("Has the presence of UNMIL increased OWN personal security?"),
covariate.labels=c("Contact with Female Peacekeeper","Contact with Male Peacekeeper","No contact with peacekeeper ","Age","Female","Savings","Muslim","High Cognitive Ability","War Trauma") ,
out="models2.html", single.row = TRUE,digits = 2, title = "Table 4: Representative Sample of Monrovia")
